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ABSTRACT 



Aims. This letter investigates the transport properties of MHD turbulence induced by the magnetorotational instability at large 
Reynolds numbers Re when the magnetic Prandtl number Pm is larger than unity. 

l\/letllods. Three MHD simulations of the magnetorotational instability (MRI) in the unstratified shearing box with zero net flux are 
presented. These simulations are performed with the code Zeus and consider the evolution of the rate of angular momentum transport 
as Re is gradually increased from 3125 to 12500 while simultaneously keeping Pm = 4. To ensure that the small scale features of the 
flow are well resolved, the resolution varies from 128 cells per disk scaleheight to 512 cells per scaleheight. The latter constitutes the 
highest resolution of an MRI turbulence simulation to date. 

Results. The rate of angular momentum transport, measured using the a parameter, depends only very weakly on the Reynolds num- 
ber: a is found to be about 7 x 10"^ with variations around this mean value bounded by 15% in all simulations. There is no systematic 
evolution with Re. For the best resolved model, the kinetic energy power spectrum tentatively displays a power-law range with an 
exponent -3/2, while the magnetic energy is found to shift to smaller and smaller scales as the magnetic Reynolds number increases. 
A couple of different diagnostics both suggest a well-defined injection length of a fraction of a scaleheight. 

Conclusions. The results presented in this letter are consistent with the MRI being able to transport angular momentum efficiently at 
large Reynolds numbers when Pm = 4 in unstratified zero net flux shearing boxes. 
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1. Introduction 

Angular momentum transport in accretion disks has been an out- 
standing issue in theoretical astrophysics for decades. To date the 
. most likely mechanism appears to be MHD turbulence driven 
' by th e magnetorotational instability (MRI, 'Balb us & HawlevI 
' [T99TI 11998). Several numerical simulations have been per- 
formed to study its properties. The most popular approach is to 
work in the local approximatio n, usin g the shearing box mo del, 
', as pioneered by |H awlev et al] (Il995h . iHawlev et al] (Il996h or 
' iBrandenburg et alJ tl995i) . These early simulations have shown 
I that MRI-powered MHD turbulence is a robust mechanism that 
transports angular momentum outw ard. The rate of transport , 
measured by the famous a parameter (Shakura & Sunvaevll973h 
depends on the field geometry but is always positive, indicating 
outward flux of angular momentum. The results obtained in the 
1990's however obviously suffered from the limited computa- 
tional resources available at that time. With no mean vertical 
magnetic field threading the she aring box (a fi eld geometry re- 
ferred to as the zero net flux casel lFromang & P apaloizou ( 2007) 
recently demonstrated with the code Zeus ( iHawlev & Stonei 
that it is indeed a problem: a decreases by a factor of two 
each time the resolution is doubled. This behavior has since been 
shown to be very robust as it has been confirmed b y simulations 
performed with codes using different algorithms (ISimon et al.l 
l2009HGuan et al.ll2009l) . This result, although it raised the con- 
cern that MRI-induced transport could vanish at infinite reso- 



lution, was interpreted as an indication that the small scale be- 
havior of the flow is an important ingredient to determine the 
rate of MRI-induced angular momentum transport: small scale 
explicit dissipation coefficients, namely viscosity and resistiv- 
ity, ne ed to be included in the s imulations. With such calcu- 
lations iLesur & Longarettil (12007 ) showed that, for a nonzero 
vertical mean magnetic field, a rises with the magnetic Prandtl 
number Pm, the ratio of viscosity over resistivity. This result 
is actually very general: it is independent of the field geom- 
etry and was also foun d for a mean toroidal magnetic field 
('Simon & Hawlev"2009^ and in the zero net flux case of interest 
here (Fromang et al. 2007.) . Recentlv ISimon et al.l (12009) mea- 
sured the numerical dissipation properties of the code Athena 
dGardiner & Stonell2008l: IStone et al.ll2008l) . They found that an 
increase in resolution amounts to an increase of the numerical 
Reynolds numbers, while keeping the eflrective magnetic Prandtl 
number (i.e. the ratio between the numerical viscosity and the 
numerical resistivity) roughly constant and equal to about two. 
In light of these results a possible interpretation of the findings of 
Fromang & Papaloizou (2007) is that a is decreasing when the 
physical Reynolds number increases at fixed Pm. If unchecked, 
this decreasing a would mean that MRI-induced MHD turbu- 
lence is ineffective at transporting angular momentum without 
a mean flux, even in systems that have Pm values higher than 
unity. Here, high resolution numerical simulations in which Re 
and Rm are simultaneously increased while keeping their ratio 
Pm constant are used to examine if this is indeed the case. 
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Model 


Resolution 


Re 




a 


i?e3125 


(128, 192, 128) 


3125 


7.9 X 10-' 


± 4.5 X 10^" 


i?e6250 


(256, 384, 256) 


6250 


5.9 X 10-^^ 


± 1.8 X IQ-'* 


^el2500 


(512, 768,512) 


12500 


8.4 X 10-^^ 


± 3.8 X 10-^ 



Table 1. Properties of the simulations and time averaged value of 
a. The errors on a are computed following Lo ngaretti & Lesuj 
(l2010l) : the time history of a is divided in N bins of size t. t is 
varied between 0.1 and 8 orbits. For each N, the standard devia- 
tion cTf^ is computed according to cr^ = {^ifli - a)IN'\^^^, where 
a; is the mean value in bin i. For large N, cr^ scales like N^^^^. 
The errors reported on a use that scaling to estimate ctn when 
T = 40 orbits, the time duration over which the mean values of 
a are calculated. 



2. Numerical setup 

In the simulations described below, the non-ideal MHD equa- 
tions (i.e. including viscosity v and resistivity ri) are solved in 
the unstratified shearing box (Goldreich & Lvnden-Bel]||1965h 
by the code Zeus ( Hawlev & Ston e 1995). The setup is identi- 
cal to that used bv iFromang et al.l (12007): the shearing box ro- 
tates around the central point mass with angular velocity Q. (thus 
defining the orbital time Tgrh = 2:7r/Q), the equation of state 
is isothermal with the sound speed cq, and the size of the box 
is fixed to {Lx,Ly,L,) - {H,nH,H), where H - co/il is the 
disk scaleheight. As mentioned in the introduction, the magnetic 
flux threading the disk vanishes in all directions. Three simu- 
lations are presented here. They share the same value for the 
magnetic Prandtl number Pm - v/rj = 4. The Reynolds num- 
ber Re = cqH/v is gradually increased from Re = 3125 (here- 
after labeled model 7?e3125) to Re = 6250 (model /;e6250) 
and finally Re - 12500 (model /?el2500). The resolution is 
increased at the same time as the Reynolds number to en- 
sure that the smallest scale features of the flow are always re- 
so lved. Model j;e3125 is identical to model 128/?e3125Pm4 
of iFromang et al.l (l2007h . for which different diagnostics have 
shown that 128 cells per scaleheight are sufficient when us- 
ing Zeus. Thus the resolutions iN„Ny,N-) = (128, 192, 128), 
(256, 384, 256) and (512, 768, 512) are adopted respectively for 
model /?e3125, /?e6250 and Rel25Q(E- For model ReUSOQ, it 
was found that early transients associated with the linear insta- 
bility kept affecting the flow for long times, resulting in pro- 
hibitively long simulations. For the computational cost of that 
simulation to remain acceptable, the following procedure was 
used: model 7?e3125 was run from f = to f = 150 orbits, start- 
ing fro m the initial state des cribed above and identical to that 
used by IFromang et al.l (120071) . At f = 60, the flow was inter- 
polated on a grid twice finer. The dissipation coefficients were 
reduced by a factor of two and the model was restarted between 
f = 60 and t - 150 orbits. This constitutes model /?e6250. This 
procedure was repeated at time f = 90 orbits to produce model 
7?el2500. The latter was run between f = 90 and t - 135 orbits. 
The properties of the three models are summarized in Table [T] 
the first column gives the label of the model, ffie second col- 
umn reports its resolution (A^i, A^,,; N,) and the third the Reynolds 
number Re for that run. All models share the same value Pm - 4. 



' Model iJe 12500, with 512 cells per scaleheight, constitutes the 
highest resolution published so far of MRI induced turbulence. With 
about 2 X 10** cells, the simulation required over 1.4 million timesteps 
to be completed and a total of about 350000 CPU hours on the CEA 
supercomputer BULL Novascale 3045 hosted in France by CCRT. 
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Fig. 1. Time history of the Maxwell stress tensor for model 
i?e3125 (blue dotted line), 7?e6250 (green dashed line) and 
7?el2500 (red solid line). The three curves are consistent with 
ffie same time averaged value for omox, independently of the 
Reynolds number. 



Finally, the last column in Table[T]gives time-averaged values of 
a that are discussed in the subsections below. 



3. Flow properties 

In the three models flow features typical of unstratified 
shearing boxes simulations are recovered: weakly non- 
axisymmetric density waves propagate radially in the box 
(Heinemann & Papaloizou 2009a b) on top of smaller scales ve- 
locity and magnetic field turbulent fluctuations, the latter exhibit- 
ing a tangled structure ty pical of Pm values higher than unity 
( Schekochihi n et a n i2004 . Below we concentrate on the trans- 
port properties of the turbulence, the shape of the kinetic and 
magnetic energy power spectra and ffie the two points correla- 
tion function. 

3. 1 . Angular momentum transport 

The angular momentum transport properties of the turbulence 
in the three models are assessed by calculating the a parame- 
ter, the sum of the Reynolds stress tensor ajiev and the Maxwell 
stre ss tensor aMax- All three coefficients are calculated as in 
Fro mang & Papaloizoul (l2007h . The time history of aMax is 
shown in Fig. [1] for models 7?e3125, 7;e6250 and ReUSOO re- 
spectively, using a dotted, a dashed and a solid line. The result is 
dramatically different from the results of lFromang & Papaloizoul 
(2007) who found without explicit dissipation a monotonic de- 
crease of OMax as the resolution was increased. Here, no such 
systematic evolution is found as Re goes up. Indeed umox ap- 
pears to vary only very weakly with the Reynolds number. This 
is confirmed by the last column of Table [T] in which the values 
of a, time-averaged between 90 and 130 orbits, are reported for 
the different models. The rate of angular momentum transport 
appears to be somewhat smaller in model /?e6250 than in model 
/?e3125 and /?e 12500. Nevertheless, the difference between the 
three simulations remains less ffian 25%. Taken together, the 
three measurements suggest that a is of the order of 7 x lO"-' 
and is fairly independent of the Reynolds number. At the very 
least, a systematic evolution of a with Re is ruled out by the 
simulations. 
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Fig. 2. Top panel: kinetic {solid line) and magnetic {dashed line) 
energy power spectrum for model /?e 12500, time averaged over 
twenty snapshots between t - 9Q and t - 120. The dotted line 
shows a power law line with index -3/2 for the purpose of 
comparison. Bottom panel: kinetic energy power spectra com- 
pensated by 1 {dotted line), 3/2 {solid line) and 5/3 {dashed 
line). Both panels are suggestive of a A;"^^^ spectrum in the range 
30 < /t < 100. 



3.2. Power spectrum 

The top panel of Fig. |2] shows the shell-averaged kinetic en- 
ergy power spectrum Ek {solid line) and magnetic energy power 
spectrum Em {dashed line) for model 7?e 12500. The latter is 
rather flat over about a decade in wavenumber (from k ~ \Q 
to k ~ 100) and is larger than the kinetic energy over that range. 
By contrast, the former displays a clear power-law behavior for 
wavenumber 20 < A; < 100. For the purpose of comparison, 
the dotted line shows a pure power-law with the index -3/2 
that nicely fits the solid line of the plot. By analogy with hydro- 
dynamic turbulence it is tempting to associate the large scale 
end of the power-law part of the spectrum with an injection 
length li„j ~ 2n/ki„i„ ~ Q.3H. Similarly, the small scale end 
can be associated with the viscous cut-off' length and is found 
to be lyisc ~ In/kyisc ~ 0.06//. This is about 32 cells at that 
resolution and is thus well-resolved by the code. Furthermore, 
results obtained in the kinematic regime of incompressible and 
homogeneous MHD turbulence suggest that the resistive length 
Ires ~ Pin^^^^lvisc (Schckochihin ct al. 2004). Thus, is of or- 
der 16 cells and also well resolved, which shows that numeri- 
cal dissipation is most likely negligible in this simulation. Given 
the still limited resolution of model /?e 12500, the reliability of 
the power-law exponent mentioned above can however be ques- 
tioned: for that purpose, the bottom panel of Fig. |2]displays three 
compensated spectra of Ek, kE^ {dotted line), k^^^Ex {solid line) 
and k^^^Ex {dashed line) respectively. First, the figure illustrates 
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Fig. 3. Top panel: plot of Ek for model /?e3125 {dotted line), 
/?e6250 {dashed line) and /?el2500 {solid line). The dotted line 
shows a power-law line with the index -3/2 for comparison. As 
Re and Rm increase, the kinetic energy display an increasing re- 
gion well-fitted by a power law with the index -3/2, while the 
viscous cut-off region moves to higher k values. Bottom panel: 
Same as the top panel, but for the quantit y kE^- On both pan- 
els the insets reproduce the results of Fromang & Pap aloizoul 
(I2007h for model STD64 {dotted line), STD128 {dashed line) 
and STD256 {solid line). 



the difficulty of a reliable determination of the exponent. Indeed, 
the power-law extends over less than a decade in wavenumber. 
Nevertheless, the dashed line, which unambiguously rises over 
the interval of 10 < k < 100, excludes a k'^^^ spectrum and 
rather suggests an exponent larger than -5/3. The dotted line on 
the other hand suggests -1 as an upper limit. Finally, the solid 
line suggests k^^^^ as a tentative fit for the power-law range of 
the spectrum (30 < k < 100). Finally, Fig. [3] (fo/j panel) com- 
pares the shape of Ek in model /?e3125 {dotted line), /?e6250 
{dashed line) and /?el2500 {solid line). For all models, the ki- 
netic energy power spectrum peaks at k ~ 10-20. For larger 
wavenumbers, the k^^^^ power-law becomes more and more ap- 
parent as the Reynolds number increases. The bottom panel of 
Fig.[3]plots the quantity kEM for the three simulations. The peak 
of each curve thus provides an estimate of the scale at which 
magnetic energy is located. It is found to lie at kpeak ~ 30- 
40, 50-60 and 70-80 respectively when Re = 3125, 6250 and 
12500. In other words, the scale at which most of the magnetic 
energy is located moves toward smaller and smaller scales as 
Rm is increased. Th is is different from the results reported by 
iHaugen et al.l (l2003l) . but not unexpected given existing theo- 
ries of small scale dynamos with large Pm (Schekochihin et al] 
2002a,b'). On both panels, the small ins ets plot the spectra ob- 
tained by iFromang & Papaloizoul (l2007h without explicit dissi- 
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Fig. 4. Structure of the correlation function ^i, in the (A^,Ay) 
plane for model 7?e3125 {left panel), 7?e6250 {middle panel) and 
7?el2500 {right panel). Its structure is only weakly dependent 
on the Reynolds number, suggesting a well-defined injection 
length. 



pation. Aside from the decrease of their amplitude with reso- 
lution, the most noticable differences with the results presented 
here are twofold: first, the kinetic energy power-spectra appear 
flatter at intermediate wavenumbers. In addition, there is more 
energy (both kinetic and magnetic) at the smallest scales of the 
box. 



3.3. Correlation length 

The shape of the kinetic energy power-spectrum described 
above suggests an injection length U„j that appears to be inde- 
pendent of Re for the range of the Reynolds numbers investi- 
gated here. However, the shell average involved in its derivation 
washes out all information about the anisotropy of the turbu- 
lence. Th is can be investigated using the two- points-correlation 
function (iGuan et al.ll20()^lDavis et al.ll2010h : 



^v{Ax) - <'Li6vi{x)6vi{x + \x)> I <l.i6v]> . 



(1) 



Here <.> stands for a volume average, 6vi{x) corresponds to the 
velocity fluctuations in the direction i and the sum is over spa- 
tial coordinates. Isocontours of in the plane Az = are shown 
in Fig. |4] From left to right, the different panels correspond to 
models Rei 125, ^ e625 and j^el2500 resp ectively. As found by 
iGuan et al.1 (l2009h and lDavis et"an (l2010l) . has an elhpsoidal 
shape for all models. The tilt angle 0,. of the major axis is 9y~% 
degrees, which agrees with the results of Guan et al. (2009). 
Following the procedure outlined by these authors (i.e. by fit- 
ting the shape of the correlation function in a given direction by 
an exponential), the correlation lengths along the major, minor 
and vertical axis of the ellipsoid are found to be (/l,„,„, A„,ax, ^;) - 
(0.08, 0.45, 0.08)//, independently of the Reynolds number. This 
suggests once more that the injection length of the turbulence is 
only weakly dependent on the Reynolds number. At the same 
time, the small difference between A„u„ or and /,,,ic as quoted 
in section [3T2I is a warning that the injection range and the dissi- 
pative range might overlap in these simulations. 

4. Conclusion 

Here zero net flux high resolution numerical simulations of 
MRI-driven MHD turbulence are used to demonstrate this re- 



sult: when Pm-A, the dependence of a on the Reynolds number 
is very weak. In all models, a ~ 1 x 10"^ to within about 15%. 
This result unambig uously shows that the de crease of a with res- 
olution reported bv lFromang & Papaloizoul (l2007h is a numeri- 
cal artifact that contains no physical information about the na- 
ture of the MHD turbulence in accretion disks. Quite differently, 
the present simulations are consistent with a nonzero value of 
a at infinite Reynolds numbers for a magnetic Prandtl number 
higher than unity. Note that this weak dependence of a with Re 
for Pm > 1 is also suggeste d by the data recently rep orted by 
ISimon & Haw lev (2009) and lLongaretti & Lesuii(l2010h respec- 
tively for a mean azimuthal and vertical magnetic field. 

In addition, a number of statistical properties of the turbu- 
lence are reported. The kinetic energy power spectrum of the 
turbulence and the two-points-correlation function of the ve- 
locity both suggest a well-defined injection length /,„y of a few 
tens of a scaleheight. For the range of the Reynolds numbers 
Re that can be probed with current resources, li„j seems to 
be independent of Re. At the highest resolution achieved here, 
the kinetic energy power spectrum displays a power-law scal- 
ing over almost a decade in wavenumber. However, given the 
limited extent of the power-law range, the precise exponent 
of this power-law cannot be accurately determined: an expo- 
nent of -3/2 appears to be consistent with the data, while a 
-5/3 exponent seems too steep. Nevertheless, as suggested in 
Sect. 13.31 the separation between the forcing and the dissipative 
scales might still be marginal. This is why a detailed compari- 
son of these exponents with exist ing MH D turbulence theories 
(IIroshnikovlll963t lKraichnan|[T96 5: Goldr eich & Sridhaii [19951) 
is probably premature at this stage. Higher resolution simula- 
tions are definitively needed. Finally, the shape of the magnetic 
energy power spectrum shows that magnetic energy is mostly 
located at small scales and shifts to smaller and smaller scales 
as Rm increases, as expect ed from small scale dynamo theory 
(Schekochihin et al. 2002a). This is consistent with the scenario 
postulated by Rincon et al. (2008) of a large scale MRI forcing 
that generates and coexists with a small scale dynamo. 
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